Skip to content

Conversation

@adcroft
Copy link
Member

@adcroft adcroft commented Dec 16, 2025

Adds CMIP diagnostics t17d and t20d.

t17d is depth of the 17 degC isotherm. 0m where the ocean is everywhere cooler that 17 degC, and sea-floor when everywhere warmer. Simialrly, t20d is the depth of the 20 degC. "depth" here is defined as distance from the sea-surface.

Result from ocean_only/gloal_ALE/z initial conditions:
image

Meta-data in output looks like:

        float t17d(Time, yh, xh) ;
                t17d:long_name = "Depth of 17 degree Celsius Isotherm" ;
                t17d:units = "m" ;
                t17d:missing_value = 1.e+20f ;
                t17d:_FillValue = 1.e+20f ;
                t17d:cell_methods = "area:mean yh:mean xh:mean time: point" ;
                t17d:cell_measures = "area: area_t" ;
                t17d:standard_name = "depth_of_isosurface_of_sea_water_potential_temperature" ;
        float t20d(Time, yh, xh) ;
                t20d:long_name = "Depth of 20 degree Celsius Isotherm" ;
                t20d:units = "m" ;
                t20d:missing_value = 1.e+20f ;
                t20d:_FillValue = 1.e+20f ;
                t20d:cell_methods = "area:mean yh:mean xh:mean time: point" ;
                t20d:cell_measures = "area: area_t" ;
                t20d:standard_name = "depth_of_isosurface_of_sea_water_potential_temperature" ;

Closes #992

@adcroft adcroft added the enhancement New feature or request label Dec 16, 2025
t17d and t20d are the depth of the 17 and 20 degC isotherms
respectively, measured from the surface. Interpolation is done
via Recon1d reconstruction and then a new function in the class
that solves that inverts the reconstruction. This solver was not
in the original class because I had originally only thought it
would be used in the new grid generator which interpolates
density and not for the reconstructed field. Since some other
diagnostics for CMIP will make use of remapping or T, I figured
we might as well use the reconstructions to do the interpolation.

- adds to new registrations and posts in MOM_diagnostics
- adds function `%x(self, k, val)` to the Recon1d class.
@adcroft adcroft marked this pull request as draft December 16, 2025 19:57
@adcroft adcroft marked this pull request as ready for review December 16, 2025 19:58
enddo
else
! Non-Boussinesq
p(0) = 0.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Setting the surface pressure to 0 will be inaccurate under ice-shelves.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point.

p(k) = p(k-1) + h(i,j,k) * ( GV%H_to_RZ * GV%g_Earth ) ! Pressure at lower interface
p(k-1) = 0.5 * ( p(k-1) + p(k) ) ! EOS fn needs pressure in the layer center
enddo
call calculate_spec_vol(tv%T(i,j,:), tv%S(i,j,:), p(0:nz-1), spv, tv%eqn_of_state)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Finding the distance between interfaces using the specific volume at the cell-center is inherently less accurate than actually using the layer-average specific volume, which is usually pre-stored in tv%SpV_avg(:,:,:). If this has not been pre-stored or there is some other reason not to use this, the layer-average specific volume can be obtained by calling average_specific_vol().

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good to know - that means we can avoid this column-wise function call. Is this available independent of USE_TEMPERATURE? How are we ensuring that tv%SpV_avg is always updated when h is updated?

xr = 1.
xo = ( t - this%ul(k) ) / slp ! First guess by regula falsi
f_at_x = this%f(k, xo)
do iter = 1,5
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please add a comment explaining why this code is using Newton's method to solve for the root of a quadratic equation rather than just using the closed-form solution. (Could it be that the profile is quadratic in T and S, but with a nonlinear equation of state the function might be nonlinear for density, for example?)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

enhancement New feature or request

Projects

None yet

Development

Successfully merging this pull request may close these issues.

CMIP7 diagnostics t17d and t20d

2 participants